/* START 3a_genHomeWork2015.do */

/* this script generates home and work locations for the 2015-2016 data for the
urban/spatial model. We adopt a before vs. after DTL2 approach:
- home: modal origin in the morning/destination in the evening
- work: modal destination in the morning/origin in the evening
- leisure: any weekend trip away from home, divided by the total number */

/* prerequisite packages: fillmissing, labutil2, gtools, parallel */
version 16.1
set more off
cap clear all
cap restore, not

/* PATHS */
do "code/config.do"

/* Initialize 8 threads */
// parallel initialize 8

/* load accessory functions */
qui do "code/99_accessoryFunctions.do"
use "make_data/2015.dta", clear

/* fix labeling from earlier */
cap label copy oa1 oa
cap label drop oa1
cap label values oa da oa

/*******************************/
/** SUBSET THE 2015-2016 DATA **/
/*******************************/

local planning_areas op dp
local subzones oa da
local units_of_analysis `subzones'
/* local units_of_analysis `subzones' `planning_areas' */
/* local units_of_analysis `planning_areas' */
local other_restrictions & (TKT_TYP_CD == 34 | TKT_TYP_CD == 38)
local suffix AdultsOnly

/* loop over all trip dates and compute home/work */
gen tripDate = dofc(START)
format tripDate %td
drop if tripDate == .

/* 0: before train line opens; 1: after */
foreach timePeriod in 0 1 {
	local qr `timePeriod'
	disp "Now computing `units_of_analysis' ///
	  flows for `qr' with suffix `suffix'"
	cap frames drop t`qr'
	/* frame put in 1/2, into(t`qr') */
	if `timePeriod' == 0 {
		frame put if (tripDate >= td("01dec2015") & ///
		  tripDate < td("27dec2015")) `other_restrictions', into(t`qr')
		}
	else {
		frame put if (tripDate > td("27dec2015")) `other_restrictions', ///
		  into(t`qr')
		}

	fv t`qr'
	if (strpos("`units_of_analysis'", "op") > 0) drop if op == . | dp == .
	if (strpos("`units_of_analysis'", "oa") > 0) drop if oa == . | da == .
	if _N > 0 {
		/* execute modal OD code */
		cap gen byte hstart = hh(START)
		cap gen byte am = hstart >= 5 & hstart <= 11
		cap gen byte pm = hstart >= 15 & hstart <= 23
		cap gen dayOfWeek = dow(dofc(START))

		/* figure out correct locations */
		cap restore, not
		cap frames drop w`qr'
		frame put if dayOfWeek >= 1 & dayOfWeek <= 5, into(w`qr')
		frame w`qr' {
			foreach timeOfDay in am pm {
				/* GENERATE MODAL HOME/WORK */
				foreach x in `units_of_analysis' {
					disp "Now generating `x'_mode_`timeOfDay'"
					cap lab drop `x'_mode_`timeOfDay'
					/* drop if `timeOfDay' == 0 */
					hashsort CRD_NUM `timeOfDay' `x'
					order CRD_NUM `timeOfDay' `x'
					parallel, by(CRD_NUM `timeOfDay'): gegen int `x'_mode_`timeOfDay' = mode(`x') if `timeOfDay' == 1, by(CRD_NUM `timeOfDay') minmode missing
					parallel, by(CRD_NUM): fillmissing `x'_mode_`timeOfDay', by(CRD_NUM)
					}

				/* LABEL VALUES */
				foreach x in `planning_areas' {
					cap label values `x'_mode_`timeOfDay' op
					}

				foreach x in `subzones' {
					cap label values `x'_mode_`timeOfDay' oa
					}
				}

				/************************/
				/* assign home and work */
				/************************/

			foreach name in pHome aHome pWork aWork {
				cap label drop `name'
				}

			/* planning areas only */
			if strpos("`units_of_analysis'", "op") > 0 {
				/* similarly specify home and work for planning areas */
				gen pHome = op_mode_am if op_mode_am == dp_mode_pm
				mdesc pHome
				
				replace pHome = op_mode_am if dp_mode_am == .
				replace pHome = dp_mode_pm if op_mode_am == . | pHome == .
				replace pHome = op_mode_am if pHome == .
				lab val pHome op
				
				gen pWork = dp_mode_am if dp_mode_am == op_mode_pm
				mdesc pWork if dp_mode_am != .
				
				replace pWork = dp_mode_am if op_mode_pm == .
				replace pWork = op_mode_pm if dp_mode_am == . | pWork == .
				replace pWork = dp_mode_am if pWork == .
				lab val pWork op

				cap frames drop work`qr'
				frame put CRD_NUM TKT_TYP_CD pHome pWork, into(work`qr')
				frame work`qr' {
					gduplicates drop
					drop if pHome == . | pWork == .
					gcollapse (count) totwork = CRD_NUM, by(pHome pWork) merge
					cap gcollapse (count) totworkw = CRD_NUM if TKT_TYP_CD == 38, ///
					  by(pHome pWork) merge
					gcollapse (mean) totwork (mean) totworkw, by(pHome pWork)
					ren pWork dp
					}
				}

			/* subzones only */
			if strpos("`units_of_analysis'", "oa") > 0 {
				gen aHome = oa_mode_am if oa_mode_am == da_mode_pm
				replace aHome = oa_mode_am if da_mode_am == .
				replace aHome = da_mode_pm if oa_mode_am == . | aHome == .
				replace aHome = oa_mode_am if aHome == .
				lab val aHome oa
				gen aWork = da_mode_am if da_mode_am == oa_mode_pm
				replace aWork = da_mode_am if oa_mode_pm == .
				replace aWork = oa_mode_pm if da_mode_am == . | aWork == .
				replace aWork = da_mode_am if aWork == .
				lab val aWork oa

				cap frames drop work`qr'
				frame put CRD_NUM TKT_TYP_CD aHome aWork, into(work`qr')
				frame work`qr' {
					gduplicates drop
					drop if aHome == . | aWork == .
					gcollapse (count) totwork = CRD_NUM, by(aHome aWork) merge
					cap gcollapse (count) totworkw = CRD_NUM if TKT_TYP_CD == 38, ///
					  by(aHome aWork) merge
					gcollapse (mean) totwork (mean) totworkw, by(aHome aWork)
					ren aWork da
					}
				}

			/* NOTE: STILL DOESN'T WORK WITH OP AND OA TOGETHER */
			if (strpos("`units_of_analysis'", "op") > 0) keep CRD_NUM TKT_TYP_CD pHome pWork
			if (strpos("`units_of_analysis'", "oa") > 0) keep CRD_NUM TKT_TYP_CD aHome aWork
			gduplicates drop
			gduplicates drop CRD_NUM, force
			}
		/* fv t`qr' */
		frlink m:1 CRD_NUM, frame(w`qr')
		if (strpos("`units_of_analysis'", "op") > 0) frget pHome pWork, from(w`qr')
		if (strpos("`units_of_analysis'", "oa") > 0) frget aHome aWork, from(w`qr')

		/* preserve */
		/* keep weekend trips only */
		if (strpos("`units_of_analysis'", "op") > 0) {
			keep if (pHome != . & pWork != .) & mod(dayOfWeek, 6) == 0
			gen leisure = pHome != dp | op == dp
			gcollapse (count) leisureTrips = dayOfWeek (mean) TKT_TYP_CD, by(CRD_NUM pHome dp)
			hashsort CRD_NUM dp leisureTrips
			}
		else if (strpos("`units_of_analysis'", "oa") > 0) {
			keep if (aHome != . & aWork != .) & mod(dayOfWeek, 6) == 0
			gen leisure = aHome != da | oa == da
			gcollapse (count) leisureTrips = dayOfWeek (mean) TKT_TYP_CD, by(CRD_NUM aHome da)
			hashsort CRD_NUM da leisureTrips
			}

		/* leisure measure 1: modal destination */
		gegen maxTrips = max(leisureTrips), by(CRD_NUM)
		gen argmax = 1 if leisureTrips == maxTrips
		/* gegen numArgmax = total(argmax), by(CRD_NUM) */

		frames put CRD_NUM TKT_TYP_CD *Home d* argmax, into(mode`qr')
		frame mode`qr' {
			keep if argmax == 1
			hashsort CRD_NUM
			by CRD_NUM: sample 1, count
			/* if repeated max, draw at random. DESTRUCTIVE STEP */
			if (strpos("`units_of_analysis'", "op") > 0) {
				gcollapse (count) totleis = CRD_NUM, by(pHome dp) merge
				cap gcollapse (count) totleisw = CRD_NUM if TKT_TYP_CD == 38, ///
				  by(pHome dp) merge
				gcollapse (mean) totleis (mean) totleisw, by(pHome dp)
				keep if dp != .
				/* fill missing planning areas */
				xtset pHome dp
				tsfill
				xtset dp pHome
				tsfill
				}

			else if (strpos("`units_of_analysis'", "oa") > 0) {
				gcollapse (count) totleis = CRD_NUM, by(aHome da) merge
				cap gcollapse (count) totleisw = CRD_NUM if TKT_TYP_CD == 38, ///
				  by(aHome da) merge
				gcollapse (mean) totleis (mean) totleisw, by(aHome da)
				keep if da != .
				/* fill missing subzones */
				xtset aHome da
				tsfill
				xtset da aHome
				tsfill
				}
			}

		/* leisure measure 2: 1/N */
		drop maxTrips argmax
		gegen totLeisCard = total(leisureTrips), by(CRD_NUM)
		gen leisFrac = leisureTrips/totLeisCard

		/* Generate leisure and workfare alt leisure counts */
		if (strpos("`units_of_analysis'", "op") > 0) {
			gcollapse (sum) totleisalt = leisFrac, by(pHome dp) merge
			cap gcollapse (sum) totleiswalt = leisFrac if TKT_TYP_CD == 38, ///
			  by(pHome dp) merge
			cap gcollapse (sum) totleiswcount = leisureTrips if TKT_TYP_CD == 38, ///
			  by(pHome dp) merge
			gcollapse (mean) totleisalt (mean) totleiswalt ///
			  (sum) totleiscount = leisureTrips (mean) totleiswcount, by(pHome dp)
			keep if dp != .
			/* fill missing planning areas */
			xtset pHome dp
			tsfill
			xtset dp pHome
			tsfill
			frlink m:1 pHome dp, frame(work`qr')
			frlink m:1 pHome dp, frame(mode`qr')
			}

		else if (strpos("`units_of_analysis'", "oa") > 0) {
			gcollapse (count) totleisalt = leisFrac, by(aHome da) merge
			cap gcollapse (count) totleiswalt = leisFrac if TKT_TYP_CD == 38, ///
			  by(aHome da) merge
			cap gcollapse (sum) totleiswcount = leisureTrips if TKT_TYP_CD == 38, ///
			  by(aHome da) merge
			gcollapse (mean) totleisalt (mean) totleiswalt ///
			  (sum) totleiscount = leisureTrips (mean) totleiswcount, by(aHome da)
			keep if da != .
			/* fill missing subzones */
			xtset aHome da
			tsfill
			xtset da aHome
			tsfill
			frlink m:1 aHome da, frame(work`qr')
			frlink m:1 aHome da, frame(mode`qr')
			}

		frget totwork totworkw, from(work`qr')
		frget totleis totleisw, from(mode`qr')
		drop work`qr' mode`qr'
		foreach var in totleis totleisw totwork totworkw ///
		  totleisalt totleiswalt totleiscount totleiswcount {
			replace `var' = 0 if `var' == .
			}
		if `timePeriod' == 0 {
			gen str period = "before"
			}
		else {
			gen str period = "after"
			}
		compress
		/* clean up: delete the non-t frames */
		cap frames drop work`qr'
		cap frames drop mode`qr'
		cap frames drop w`qr'
		if (strpos("`units_of_analysis'", "op") > 0) ///
		  save "make_data/`qr'`suffix'.dta", replace
		if (strpos("`units_of_analysis'", "oa") > 0) ///
		  save "make_data/a`qr'`suffix'.dta", replace
		} /* _N > 0 check ends here */

	frame change default
	}

/* APPEND 2015-2016 DATA TO EXISTING PANEL DATA FROM `3_genHomeWork.do` */

cap frames drop collated
fv collated
if (strpos("`units_of_analysis'", "op") > 0) {
	use "make_data/panelODByWeek`suffix'.dta", clear
	foreach timePeriod in 0 1 {
		cap append using "make_data/`timePeriod'`suffix'.dta"
		}
	save, replace
	}

if (strpos("`units_of_analysis'", "oa") > 0) {
	use "make_data/panelODAByWeek`suffix'.dta", clear
	foreach timePeriod in 0 1 {
		cap append using "make_data/a`timePeriod'`suffix'.dta"
		}
	save, replace
	}
